Correlated electron-phonon transport from molecular dynamics with quantum baths 
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Based on generalized quantum Langevin equations for the tight-binding wave function amplitudes 
and lattice displacements, electron and phonon quantum transport are obtained exactly using molec- 
ular dynamics (MD) in the ballistic regime. The electron-phonon interactions can be handled with 
a quasi-classical approximation. Both charge and energy transport and their interplay can be stud- 
ied. We compare the MD results with that of a fully quantum mechanical nonequilibrium Green's 
function (NEGF) approach for the electron currents. We find a ballistic to diffusive transition of 
the electron conduction in one dimensional chains as the chain length increases. 

PACS numbers: 05.60.Gg, 72.10.Bg, 63.20.kd, 73.63.-b 
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I. INTRODUCTION 

The interaction of electrons with phonons in open 
nonequilibrium molecular structures is of great impor- 
tance within the context of molecular electronics^ 2 -. A 
variety of methods at different levels of sophistication has 
been used to study this problem, each working at a spe- 
cific parameter ranged. The perturbative approach with 
a self-consistent Born approximation (SCBA) works well 
when the electron-phonon interaction (EPI) is weak, and 
has been used in the first-principles study^. In the strong 
interaction limit, it is possible to eliminate the bilinear 
EPI term via a canonical transformation 5 . This latter 
approach has only limited use in a minimum model cal- 
culation, where there is only one single electron degree of 
freedom (DOF) interacting with one single phonon DOF. 
It is also possible to study the coherent electron-phonon 
dynamics in the full coupling regimes using the scattering 
theory^, but this kind of methods ignores dephasing be- 
tween electrons and phonons. Hybrid approaches exist, 
where the electron part is treated quantum-mechanically, 
while the phonon system is handled by classical MD 7 - 
with quantum corrections^. Most of the above methods 
are developed within the context of electronic transport. 
The inclusion of phonon transport appears only very re- 
cently, mainly using the NEGF approach*^. 

Molecular dynamics is usually viewed as a method 
that produces only classical results. In this paper, we 
introduce a new MD method to study the correlated 
electron and phonon transport in open molecular junc- 
tions for the quantum systems. It is based on a gener- 
alized Langevin equation 11 for electrons and phonons, 
which so far have been used to study their quantum 
transport separatel y 12 ' 13 . The formalism is exact in 
the ballistic case, i.e., without the EPI. Quasi-classical 
approximation^ is made to the full quantum many-body 
problem for interacting systems. It does not have to as- 
sume a bilinear form of the EPI Hamiltonian, and it is ap- 
plicable to the full electron-phonon coupling range. More 
importantly, the method can simulate large systems. In 
the rest of the paper, we introduce a model system, de- 
rive the quantum Langevin equations, and analyze the 



approximation involved. We present the MD numerical 
results of molecular chains, and compare with those from 
NEGF method. 



II. MODEL AND THEORY 

Consider a typical LCR structure for transport study, 
where a molecular structure (C) is connected with two 
semi-infinite leads (L and R) as electron and phonon 
reservoirs. The two leads are linear systems in their re- 
spective thermal equilibrium states characterized by the 
chemical potential and temperature. Possible manybody 
interactions only exist in the central region. The total 
Hamiltonian is the sum of the two subsystems and their 
interaction, H c + i? p h + H cpl . The phonon part is 

H ph = £ H« h Mu L ) T v£u c ^u c ) T V™u R +V n , (1) 



a=L,C,R 



where H^ h 



\{u a ) T K a i 



is a col- 



umn vector consisting of all the displacement operators 
in the a region, and u a is its conjugate momentum. The 
atomic mass has been absorbed into m — xffru 



K a 



is the spring constant matrix. is the coupling ma- 

trix between the left lead and the central molecule, and 
— (V^ C ) T , similarly for V£ R . V n is an anharmonic 
potential, which only depends on vF . The electron sub- 
system is given in a tight-binding form in an orthogonal 
basis, 



a=L,C,R 



E 

a=L,R 



c c W n Ca c a 



h.c 



(2) 



c a (c al ) is the column (row) vector containing all the an- 
nihilation (creation) operators in the a region. V^ a has 

a similar meaning as V^ Q , and V^ a = {V" c )^ ■ h.c. rep- 
resents Hermitian conjugate. The total electron energy 
under the Born-Oppenheimer approximation depends on 
the position of the atoms, so that we can make a Taylor 
expansion of it about the atomic equilibrium positions, 



2 



and obtain the electron-phonon interaction terms (e.g., 
from a first-principles calculation) 



and 



dM^c.UkUi + 



ijk 



2 ^ 

i,j,k,l 



L %3 ^3 



(3) 



H cp i includes all the higher order terms of the Taylor 
expansion. The superscript C has been omitted since 
EPI only takes place in the center part. M*j and MM are 
the first and second order EPI coefficients, respectively. 

Working in the Heisenberg picture, we obtain the equa- 
tions of motion for operators u a and c Q , e.g., for c, 



T a c a + V c aC c c , (a = L, R), 



i c 

it c = T c c c + VF L c L + V B CR c R + \c c ,H t 



(4) 
(5) 



We set h = 1, e = 1 throughout the formulas. The lead 
operators can be solved formally, 

J*{t)=ig r a (tM)c a {tx)+ f 9 r a (t,t')V e aC c c (t')dt' : (6) 

where g r a {t,t') = -i9(t - t')([c a (t), J a {t')] + ) is the elec- 
tron retarded Green's function for the lead a. It satisfies 

^9 r a (t, + g r a (t, t)T° = -IS(t t% (7) 

with the boundary condition g^(t, t') = (t < t'). Using 
Eq. ©, the equation of motion of the central operator 
reads 



ic 



c 



T c c c + f Yr{t,t')c c (t l )dt'+£+y"M k u k c c 



(8) 



Similar equation can be derived for the phonon displace- 
ment operators^, 

u c = -K c u c + F n - J t \ TF(t, t')u c (t')dt' + j] 



M& 



(9) 



F n is the force due to anharmonic effect. The last terms 
of Eqs. (HI) and (JU) are due to EPI. We have only kept 
the first order term of the Taylor expansion, although 
inclusion of higher orders is straightforward. Equations 
© and © have the form of the generalized Langevin 
equation for the quantum Brownian motion^. 

Let us try to understand these two equations. The 
damping kernels S r = E£ + YT R and IT = W L + IIJj 
are the electron and phonon retarded self-energies in the 
NEGF formalism. They are defined as, e.g., for electron 



X r a (t,ti) = V c Ca g r a (t,ti)V c aC , (a = L,R). (10) 

In the wide-band limit, the coupling with the leads 
does not depend on the energy. The damping kernel 
approaches memoryless <5-function in the time domain. 
£ = £i(i) + £r(£) and 77 = r/i(t) + ilR.{t) are electron and 
phonon random noises due to the leads (a = L, R) 



Ut)=^ a g r a (t,t 1 )c a (h), 



(11) 



Va(t) = V^ a [d r a (t,ti)u a (h) - <£(MiK(*i)J ■ (12) 

d r a (t,tx) = -i6(t - ii)([u Q (t),u Q (^i) T ]) is the lead re- 
tarded Green's function for phonons. In the leads the 
electron and phonon subsystems do not couple. They are 
both linear systems. In addition, the left and right lead 
are completely independent. The statistical properties 
of the random noises are determined by the equilibrium 
ensembles at the remote pass, t\. Working in the eigen- 
mode representation, we can show that the expectation 
value of each noise term is zero. We can also obtain their 
correlation matrices, e.g., for electrons 



(13) 



As expected, it does not depend on the initial time t\, 
and is time translationally invariant. It is convenient to 
work in the Fourier domain, 



E a (t-t')e Mt - t ' ) dt = f°(uj)r°[uj}. (14) 



fe(<jj) is the Fermi distribution function. r™[tj] = 
i(E^[w] — S°[w]) denotes the coupling with the leads. 

is positive semi-definite, as required from a clas- 
sical noise correlation. The phonon noise has a similar 
relation. A symmetric form is used herein 



+ 00 



F a \uj] 



(15) 



= ( /pi» + 2 ) r phM 



where /^(lj) is the Bose distribution for phonons, T p h is 
similar to T c . 

We notice that the noise Eqs. (TTTj) and (|12[) contain op- 
erators, which satisfy anti-commutation or commutation 
relations. Electrons and phonons need different treat- 
ment. Equations (|13|) and (TT4")) are only applicable to 
electrons. To study the hole transport, we need to use 
the correlation matrix (£a W£tU^'))- For phonons a sym- 
metrization is needed to eliminate an imaginary part of 
the correlation. In both cases the relation between the 
damping and the noise term is a kind of manifestation of 
the quantum fluctuation-dissipation theorem. 

The electrical and energy current can be obtained from 
different methods. We can use the current continuity 
condition. In the case of a discrete Hamiltonian, the 
electrical current from cell j — 1 to cell j is, with only the 
lowest EPI term included, 

'('•//,, :<■.! : -X'V U L I''' "'^ h - ( '-)- (16) 

k 

We can also get the current from each lead by studying 
the time derivative of the electron number 



dN a 
dt 



i(c cj/ B a - h.c), (17) 
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where B a = V Ca c a =£ a + j£ E r a (t,t')c c (t')dt' . In the 
same way, the electron energy current is 



IS = -^ = -(B^ c +h.c). 



(18) 



So far the formal quantum Langevin equations are in 
terms of operators. To perform a MD simulation, we need 
to turn the operators into numbers. This is achieved by 
taking their quantum mechanical expectation values at 
the beginning of the dynamics. It is reasonable to assume 
that the central region and the two leads are decoupled at 
that time. The two baths assume canonical equilibrium 
distributions, and the central region is in an arbitrary 
state denoted by the density matrix p c . The expectation 
value of any operator A c is (A c ) = Tr{ ( o c A c }. Tak- 
ing the expectation value of these operators, generating 
the noise series using their correlations^, the operator 
Langevin equations are turned into c- number equations. 
For products of operators, mean-field type approxima- 
tion is used, e.g., (cu) w (c)(u). MD simulation can be 
done using these two equations. The final result is the 
ensemble average over the initial states. To evaluate the 
electrical current, the operators in Eqs. (fTfJl - fLS]) are re- 
placed by the c-numbers got from MD simulation, and 
also c' replaced by c*, which is the complex conjugate of 
c. By doing this, we have taken the classical approxima- 
tion to the operators. 

One may cast doubt that this approximation may be 
too inaccurate to give reasonable results for the fermionic 
system for the electrical current. However, we can 
show rigorously that for the ballistic case the classical 
Langevin dynamics with the appropriate noises gives ex- 
actly the same result as that predicted by the NEGF 
metho d 13 ^ 17 ' 18 . To do this, we write the Langevin equa- 
tions in the frequency domain 



c c \cu] 



u c \cu] 



= G) 




+ J M k u k luj>]c c luj-uj>]^pj(l9) 



F„ 



■ c \u;]Mc c [u;-aj']^pj. 



We also have 



(20) 



(21) 



In the ballistic case, we write Eq. (fT7|) in the energy 
domain and substitute Eqs. (|19H2"Tj) into it. Using the 
noise correlation Eq. (|14p . after some rearrangement, we 
get exactly the Meir-Wingreen formula in the NEGF 
method. In the presence of EPI, Eqs. (|H?1 - |2"U|) are cou- 
pled. Repeated iteration with respective to c c [oj] and 
u c [Li;] gives an infinite series of terms. Analysis of these 
terms shows that the quasi-classical approximation in- 
cludes both the crossed and the non-crossed Feynman 
diagrams. But it only reproduces correctly part of the 



high order terms in the series, e.g., out of the seven lower- 
est order nonlinear self-energy graphs, two of the graphs 
involving G > is replaced by — G < . These wrong terms 
are not important when the electron number per site in 
the center region is small or the EPI is not strong, which 
defines the application range of the quasi-classical ap- 
proximation. 



III. NUMERICAL RESULTS AND 
DISCUSSIONS 

To illustrate the present approach, we take a simple 
one-dimensional (ID) atomic chain connected with two 
ID leads and simulate the coupled equations {8]| and ([9]) 
on computer. Each atom has only one displacement de- 
gree of freedom and one spinless electron state. We take 
the two leads to be the same with spring constant hi, hop- 
ping matrix element —hi, and electron onsite energy E\. 
k c , —h c , and e c denote those of the central part. Their 
couplings are — v e and — v p h for electrons and phonons. 
Some of the matrices, e.g., T c and V e r , are given by 




V LC = 




(22) 



(23) 







The lead Green's functions have analytical solutions 10 . 
The anharmonic force F n is turned off in order to perform 
a comparison with the NEGF method. The voltage is 
applied by shifting the chemical potentials of the two 
leads. A tight-binding SSH type EPI terrr^ 



H, 



epi 



L-l 
i=l 



(cjcj+i + ct +1 Cj) (u i+ i - Ui) (24) 



is used in the simulation. The Langevin equations, with 
all the operators replaced by their expectation values, 
are numerically solved using a fourth order Runge-Kutta 
method. A time-step of At = 5 x 10~ 17 s and 10 6 
MD steps are used for each data point. As for the 
NEGF results, the Meir-Wingreen expression for electri- 
cal current 18 , I a = £ /Tr{G>S< - G<E>}dw, is used. 
The greater (lesser) self-energy S>[w] (£<[w]) is due to 
the lead a. G > [u>] (G < [oj\) is the greater (lesser) Green's 
functions of the central region. A finite difference is used 
to calculate the quantum conductance from the electrical 
current. 

We first demonstrate that the MD and the NEGF 
method give the same results in the ballistic case. Fig- 
ure [T] shows the ballistic electron conductance of a two- 
atom chain as a function of the hopping matrix element 
between them h c . When h c = 0.1 eV, the conductance 
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FIG. 1: Ballistic electron quantum conductance as a function 
of hopping matrix element between the two atoms ft c at 1 K. 
Other parameters are e c = ei = 0, hi =0.1 eV, v e =0.1 eV. 
The line is from NEGF and the dots are MD. 



reaches a maximum value corresponding to one quantum 
unit (e 2 /2ttK). Going apart from this value in both di- 
rections leads to conductance decrease. The MD and the 
NEGF method give exactly the same results within the 
statistical errors of the MD simulation. This can also be 
seen from the ballistic I-V curve in Fig. (the upper 
curve in the main panel). 

Now we turn on the EPI. The main panel of Fig. [2] 
shows the I-V characteristics of the two-atom junction. 
The lower and upper curves are with and without EPI, 
respectively. In the presence of EPI, both methods give 
approximate results. The NEGF results are based on the 
SCBA, where only the non-crossed Feynman diagrams 
are included in the self-energies^. The MD method 
is non-perturbative and includes both crossed and non- 
crossed diagrams, but only part of these diagrams are 
treated correctly. As a result, the electrical current from 
MD is lower than that from the NEGF method. This can 
also be seen from the inset of Fig. [2] where we change the 
electron-phonon interaction streghth at an applied bias 
of 0.2 V. 

A detailed analysis of the high order terms in the quasi- 
classical approximation shows that its accuracy depends 
much on the electron average occupation number in the 
center region. When the electron number in small, the 
diagrams that the quasi-classical approximation treats 
incorrectly is not important. In this regime, the MD 
method should be accurate quantitatively. Out of this 
regime, it can only gives qualitatively results. These anal- 
ysis is confirmed in Fig. [3J where we show the electrical 
current and average electron number per atom as a func- 
tion of the electron onsite energy in the center region. 
The electron number from the two methods shows slight 
discrepancy only when the onsite energy is very low. The 
MD electrical current agrees with the NEGF method only 
when the electron number is below 0.3. 

The MD approach has its advantage: it can handle 
much larger systems than the NEGF method. This is 



FIG. 2: Current-voltage characteristics of the two-atom chain 
at 1 K with the following parameters: hi = 1.0 eV, h c = 0.1 
eV, v e = 0.32 eV, e c = £; = 0, ki — k c — 0.5 eV/(amu 
A 2 ), v ph = 0.1 eV/(amu A 2 ), and m = 0.2 eV/(amu5A). 
A small onsite spring constant fco = 0.2k c is applied for the 
whole structure. MD results are shown in points, and NEGF 
in lines. The filled dots and the straight line are the ballistic 
results. The lower line and the unfilled dots are results with 
EPI. The inset shows the electrical current as a function of 
EPI strength m at V = 0.2 V. 
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FIG. 3: The electrical current I (a) and the average electron 
number per site n (b) at V = 0.04 V as a function of electron 
onsite energy e c in the center. The electron-phonon interac- 
tion m = 0.2 eV/(amu2 A). All other parameters are the same 
with Fig. H 

easy to understand. Given the total degrees of freedom 
N, we only need to solve a set of 2N coupled equations 
in the MD method. While in the NEGF method matrix 
multiplication and inverse need much longer computer 
time (of order TV 3 ). In Fig. HJ we show the length de- 
pendence of the electron conductance. Study of this ef- 
fect using the NEGF method is formidable. Besides the 
long computer time needed, convergence is also hard to 
achieve for long chains. From the log scale plot, we find 
a length independent conductance for short chains and 
close to inverse linear (1/L) dependence for long chains. 
This corresponds to a ballistic to diffusive transition of 
the electronic transport. This transition takes place ear- 
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FIG. 4: Log scale plot of the electron conductance as a func- 
tion of chain length for m — 0.05 eV/(amu2A), hi = h c — 
v e = 0.1 eV, and e c — £i — 0. Phonon parameters are the 
same with Fig. [2] 



lier at 300 K due to more available phonons for scattering. 
Previous study of this transition relies on a phenomeno- 
logical method- 2 -. Thus a first-principle method that is 
able to cover both regions is highly desirable. The MD 
method proposed here could be one candidate. 

The NEGF results should be valid at small values of 
EPI. But at intermediate interaction range, no good ap- 
proximation exists. From this point of view, the MD 
method proposed here provides an alternative nonper- 
turbative way to study the correlated electron-phonon 
dynamics in the intermediate EPI regime, although it 
is only quantitatively accurate when the electron oc- 
cupation number is small. The MD method does not 
depend on the forms of the EPI Hamiltonian and the 
phonon anharmonic potential, though not exploited here. 
More importantly, it can handle much larger systems 
than the NEGF method. Further improvement of the 



results may be obtained by including higher order quan- 
tum corrections 8,21 . 



IV. CONCLUSIONS 

In summary, we have proposed a MD method to study 
the correlated electron and phonon transport in open 
nonequilibrium molecular structures. It is based on the 
generalized quantum Langcvin equations. The effects 
of the leads are reflected in the Langevin equations as 
noises and damping terms, which satisfy the quantum 
fluctuation-dissipation theorem. Quantum effects of the 
leads are taken into account properly at least for the 
electrical or energy current calculation. The method 
gives exact results for both electrons and phonons in 
the ballistic transport regime. When there is EPI, it 
is a quasi-classical approximation. The approximation is 
valid when the electron occupation number in the cen- 
ter region is small. The method shows its advantages in 
treating large systems, where fully quantum-mechanical 
study is formidable. We illustrate this by studying the 
ballistic to diffusive transition of the electrical conduc- 
tance in ID chains. Although only examples of electrical 
currents are presented here, it has other applications. For 
example, we can also study thermoelectric transport in 
molecular structures. 



Acknowledgments 

The authors thank Per Hedegard, Mads Brandbyge, 
Jian Wang, Lifa Zhang, and Yong Xu for discussions. 
This work was supported in part by a Faculty Research 
Grant (R-144-000-173-101/112) of National University of 
Singapore. 



* Electronic address: tower.lu@gmail.com 

1 M. Galperin, M. A. Ratner, and A. Nitzan, J. Phys.: Con- 
dens. Matter 19, 103201 (2007). 

2 N. Agrait, A. L. Yeyati, and J. M. van Ruitenbeek, Phys. 
Rep. 377, 81 (2003). 

3 A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 
245302 (2004). 

4 T. Frederiksen, M. Brandbyge, N. Lorente, and A. -P. 
Jauho, Phys. Rev. Lett. 93, 256601 (2004). 

5 G. D. Mahan, Many-Particle Physics (Kluwer Aca- 
demic/Plenum Publishers, 2000), 3rd ed. 

6 H. Ness and A. J. Fisher, Phys. Rev. Lett. 83, 452 (1999). 

7 C. Verdozzi, G. Stefanucci, and C.-O. Almbladh, Phys. 
Rev. Lett. 97, 046603 (2006). 

8 A. P. Horsfield, D. R. Bowler, A. J. Fisher, T. N. Todorov, 
and C. G. Sanchez, J. Phys.: Condens. Matter 17, 4793 
(2005). 

9 M. Galperin, A. Nitzan, and M. A. Ratner, Phys. Rev. B 
75, 155312 (2007). 

J. T. Lii and J.-S. Wang, Phys. Rev. B 76, 165418 (2007). 



11 G. W. Ford, M. Kac, and P. Mazur, J. Math. Phys. 6, 504 
(1965). 

12 D. Segal, A. Nitzan, and P. Hanggi, J. Chem. Phys. 119, 
6840 (2003). 

13 A. Dhar and B. S. Shastry, Phys. Rev. B 67, 195405 (2003). 

14 A. Schmid, J. Low Temp. Phys. 49, 609 (1982). 

15 J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007). 

16 P. Hanggi and G.-L. Ingold, Chaos 15, 026105 (2005). 

17 C. Caroli, R. Combescot, P. Nozieres, and D. Saint- James, 
J. Phys. C : Solid State Phys. 4, 916 (1971). 

18 A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 
50, 5528 (1994). 

19 W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. 
Lett. 42, 1698 (1979). 

20 S. Datta, Electronic Transport in Mesoscopic Systems 
(Cambridge University Press, 1997). 

21 O. V. Prezhdo and Y. V. Pereverzev, J. Chem. Phys. 113, 
6557 (2000). 



